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ABSTRACT 
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This  paper  deals  with  fitting  two-dimensional  stationary 
random  field  (RF)  models  to  images.  We  assume  that  the  given 
image  is  represented  on  a  torus  lattice,  obeyinq  an  R.F.  model 
driven  oy  uncorrelated  noise.  The  stochastic  model  is  cnaracterizea 
by  a  set  of  unknown  parameters.  We  describe  two  sets  of  exper¬ 
imental  results.  First,  by  assigning  values  to  parameters  in 
the  stationary  range,  two-dimensional  patterns  are  generated. 

It  appears  that  quite  a  variety  of  patterns  can  be  generated. 

Next  we  consider  the  problem  of  estimating  the  parameters, 
given  an  arbitrary  image.  By  assuming  a  Gaussian  structure  for 
the  noise,  we  given  an  iterative  scheme  to  estimate  the  unknown 
parameters.  We  also  implement  a  decision  rule  to  choose  an 
appropriate  set  of  neighbors  for  the  image.  The  theory  is 
illustrated  by  applying  it  to  synthetic  patterns. 
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1.  Introduction 

Random  field  models  have  many  applications  in  image  pro¬ 
cessing  and  analysis;  for  instance,  they  can  be  used  for  the 
design  of  image  enhancement  or  restoration  algorithms  [1-3]  , 
for  image  coding  [4-6] ,  and  for  characterization  of  textures 
[7-8].  Typically,  an  image  is  represented  by  two-dimensional 
scalar  data,  the  gray  level  variations  defined  over  a  square 
grid.  One  of  the  important  characteristics  of  such  data  is  the 
statistical  dependence  of  the  gray  levels  within  a  neighbor 
set.  For  example,  y(i,j),  the  scalar  gray  level  at  position 
(i,j),  might  be  statistically  dependent  on  the  gray  levels  over 
a  neighbor  set  that  includes  { (i-1, j) , (i+1, j) , (i, j-1) , (i, j+1) } . 
This  is  in  contrast  to  the  familiar  time  series  models  where 
the  dependence  is  strictly  on  the  past  observations. 

Image  models  which  include  dependence  in  all  directions 
(referred  to  as  neighbor  set  models  in  the  sequel)  have  been 
considered  recently  [3,6, 7, 8].  The  neighbor  set  dependence 
might  include  the  four  nearest  pixels  (east,  west,  north,  and 
south  neighbors),  the  eight  nearest  pixels  [3,8],  or  all  the 
pixels  inside  a  square  window  surrounding  the  pixel  at  (i,j) 

[6].  In  these  models,  the  observation  y(i,j)  is  written  as 
a  linear  weighted  sum  of  the  observations  over  the  corresponding 
neighbor  set  and  an  uncorrelated  noise  sequence  and  is  charac¬ 
terized  by  a  set  of  coefficients  and  the  variance  of  the  noise 
driving  the  model. 


Prior  to  the  use  of  these  models,  two  problems  have  to  be 
tackled,  namely  the  estimation  of  the  unknown  parameters  and 
the  choice  of  an  appropriate  neighbor  set  for  the  given  image. 

The  second  problem  has  received  some  attention  in  the  litera¬ 
ture  [8-11] .  The  parameter  estimation  is  usually  handled  by 
the  maximum  likelihood  (ML)  method.  This  involves  imposing 
a  Gaussian  structure  on  the  noise  and  deriving  an  expression 
for  the  likelihood  of  the  observations.  Unlike  the  cases  of 
one-dimensional  time  series  models  or  two-dimensional  causal 
models,  deriving  an  expression  for  the  likelihood  of  the  obser¬ 
vations  poses  some  difficulties  for  RF  models.  This  is  essen¬ 
tially  due  to  the  fact  that  the  Jacobian  of  the  transformation 
matrix  from  the  noisy  variates  to  the  observations  is  not  unity 
and  is  difficult  to  evaluate.  Whittle  [12]  developed  an  asymp¬ 
totic  approximation  for  the  determinant  term,  and  developed  approxi 
mate  expressions  for  the  likelihood  function.  Using  spectral 
representation  of  the  RF,  likelihood  functions  in  the  transform 
domain  were  considered  in  [9,11].  The  problem  of  evaluating 
the  determinant  can  be  avoided  by  making  assumptions  about  the 
representation  of  the  lattice.  Specifically,  by  assuming  repre¬ 
sentation  on  a  torus  lattice  [8,10]  (where  the  image  is  assumed 
to  be  folded  over  a  torus) ,  explicit  expressions  can  be  derived 
for  the  determinant  term,  as  the  transformation  matrix  possesses 
a  block  circulant  structure  whose  eigenvalues  can  be  written  down 
explicitly.  We  use  this  torus  representation  developed  in  [8,10] 


and  write  explicit  expressions  for  the  likelihood  of  the 
observations.  Since  the  likelihood  function  is  nonquadratic 
in  the  parameters,  the  estimates  have  to  be  determined  by  using 
numerical  optimization  procedures  such  as  Newton- Raphson,  etc. 

[13] .  To  save  some  computational  effort,  we  suggest  an  itera¬ 
tive  scheme,  defined  by  using  a  logarithmic  approximation  to 
the  determinant  term.  This  method  yields  estimates  that  are 
close  to  ML  estimates. 

The  second  problem  considered  in  fitting  RF  models  is  the 
choice  of  appropriate  neighbors  in  images..  From  one-dimensional 
time  series  analysis  it  is  known  that  the  use  of  an  appropriate 
model  leads  to  good  results  in  forecasting  and  similar  applica¬ 
tions.  The  problem  of  choice  of  appropriate  neighbors  has  been 
considered  in  the  literature  [9,11].  The  derivation  of  asymp¬ 
totically  consistent  decision  rules  for  the  above  problem  is 
given  in  [11] ;  they  are  based  on  the  corresponding  decision 
rules  for  discriminating  between  different  autoregressive  models 

[14] .  We  implement  this  decision  rule  for  choosing  between 
different  neighbor  sets. 

The  usefulness  of  the  estimation  scheme  and  the  decision 
rule  for  the  choice  of  neighbors  is  demonstrated  by  applying 
them  to  synthetic  patterns,  the  underlying  true  model  of  the 
synthetic  pattern  being  known.  This  leads  us  to  the  problem 
of  generating  synthetic  patterns.  Computationally  elegant 


solutions  using  torus  representations  for  generating  synthetic 
patterns  have  been  developed  in  [10] .  We  use  this  scheme  and 
generate  two-dimensional  patterns.  The  patterns  are  quite 
varied  and  display  the  role  played  by  the  neighbor  sets  consi¬ 
dered  and  the  values  of  the  coefficients. 


The  organization  of  the  paper  is  as  follows:  In  Section  2, 
we  consider  the  estimation  problem  and  develop  an  iterative 
method.  The  implementation  of  the  decision  rule  for  the  choice 
of  appropriate  neighbors  is  also  discussed.  In  Section  3,  we 
give  the . experimental  results.  Discussion  is  presented  in 
Section  4. 
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2.  Estimation  of  parameters  in  RF  models 

Assume  that  the  observations  y(i,j),  (i,j)fcft,  have  zero 
mean  (which  can  be  usually  achieved  by  subtracting  the  sample 
mean  from  the  original  observations)  and  that  y ( - )  obeys  the  RF 
model  in  (2.1),  the  neighbor  set  being  denoted  by  N: 

y(i,j)  =  E  <>y  ( (i,  j)  +  (k,£) )  +/p<*>(i,j),  (i,j)tft  (2.1) 

(k ,1)  €N  K' 

In  (2.1),  (0,p)  are  unknown  parameters  and  <*>(•)  is  an  i.i.d. 
noise  sequence  with  zero  mean  and  unit  variance.  Typically, 
the  neighbor  set  N  could  include  dependence  on  nearest  neighbors 
in  the  north,  east,  south,  and  west  directions,  denoted  as 
{  (-1,0) , (0,1) , (1,0) , (0,-1) } .  To  ensure  stationarity ,  the 
coefficients  (6^  £,  (k,£)€N>  must  obey 

|  E  9V  »z*  z^|  <1  whenever  |z. |  =  [z,l  =  1  (2.2) 

(k ,1)  fcN  K,c  1  z  '  1  * 

Since  we  are  working  with  a  finite  image,  for  an  arbitrary 
neighbor  set  N,  the  neighbors  of  boundary  pixels  are  not  defined. 
Hence  the  image  is  assumed  to  be  folded  into  a  torus  so  that 
(2.3)  is  valid  for  all  (i,j)e«: 

y  t  (i,  j)  +  (ix,  jx)  ]  *  yl  (i+ij-IJmod  M+l,  (j+j^-Dmod  M+l]  (2.3) 
The  torus  assumption  ensures  that  all  the  relevant  neighbors  of 
any  y(s)  belonging  to  the  finite  image  are  well  defined. 

Letting  yT  *  (y  (1,1) , . . .  ,y  (1,M) , . . .  ,y  (M,M) )  ,  wT  =  (<*>(1,1),..., 
<*)  (1,M)  ,...,<*>  (M,M) ) ,  (2.1)  can  be  rewritten  as 


where  B(0)  is  a  block  circulant  transformation  matrix  from  the 
noisy  variates  to  the  observations  with  a  typical  structure 
such  as 


B  (0)  = 


i,i 

B1 , 2  *  *  * 

B1,M 

1  ,M 

B1  1  ... 

bi,m-i 

.  Bl#l 


(2.5) 


For  instance  when  the  neighbor  set  of  dependence  N  is 


{(-1,0), 

(0, 

rl),  d,0). 

(0,-1)} 

we  have 

Bl.l 

= 

circulant 

(i,-0Oi 

,1,°' •• • '"60, 

Bl,  2 

= 

circulant 

<-6i,o< 

r0,...,0) 

bi,m 

= 

circulant 

(— q»0,...,0) 

and 

= 

0  j  * 

1,2,M 

For  notational  convenience  the  sum  over  (i,j)fcft  is  simply 
denoted  by  ft. 

Given  an  image,  our  interest  is  to  estimate  the  parameters 


characterizing  the  image.  One  of  the  popular  methods  of  estima¬ 
tion  is  the  classical  least  squares  method.  This  method  yields 
the  estimates  defined  below: 


(2.6) 


$  =  Arg  min  {J.  (0) } 

OfcM  1  ~ 

where 

J, (6)  =  E(y(i,j)  -  0Tz(i,j))2  (2.7) 

T  0 

and  z  (i,j)  =  (y [ (i, j) + (k,f) 3 ,  (k,£)fcN} 

By  performing  the  minimization  in  (2.6)  we  have 

$-  [£z(i#j)zT(i,j)]-1(Zz(i,j)y(i,j))  (2.8) 

a"  ~ 

and 

p  -  -4  Z(y(i,j)-0Tz(i, j))2  (2.9) 

M  n 

One  of  the  drawbacks  of  this  method  is  that  in  general  $ 
is  not  always  consistent  [12,15].  Another  popular  method  of 
estimation  is  the  maximum  likelihood  (ML)  method  which  yields 
asymptotically  consistent  and  efficient  estimates.  To  obtain 
an  expression  for  the  log  likelihood  function,  we  impose  a 
Gaussian  structure  on  the  noise  sequence  <*>(•).  Using  the 
Gaussian  assumption,  the  likelihood  of  the  observations  can 
be  written  as 

*np(y|6,p)  =  In  det  B(0)-(M2/2)fn  2irp  -  ±  E  (y  (i ,  j ) -0Tz  (i , j  ) )  2 

a 

(2.10) 

From  [10],  due  to  the  block  circulant  structure  of  B(0),  we 
have 

det  b(0)  -  nu-eTi.)  (2.ii) 

n  ~ 

where 

.  .  k(i-l)+f (j-l) 

!ij  1A0 


,  (k,f)6N  } 


and 


XQ  *  expt/^^n/M} 

Using  (2.11),  we  get 

In  p(y j  8,p)  *  Z£n(l-0Ty.  .)-(M2/2)£n  2irp  -  Z (y (i , j ) -eTz (i, j ) ) 2 


-  ~i3 


2p  n 


(2.12) 


To  avoid  computations  involving  complex  quantities  in  (2.12), 
we  use  the  following  lemma. 


Lemma  1: 


ft 

where 


Ztn  (1— 0 A 4* .  . )  =  0.5  Zfn (1-28TC . .+0TQ.  .6) 


~  ~X3 


-  -  ~x]~ 


(2.13) 


?ij  =  Vcos(X0((i-l')k+(j-lU)],‘  (k,f)tN} 
S*.  =  (sin[X0((i-l)k+(j-l)t)],  (k,f)tN> 


Q.  .  =  C.  .cT.  +  S.  .sT.  (2.14) 

~lj  -ID-13  -13-13  v  ' 

Using  Lemma  1,  the  likelihood  function  can  be  written  as 

In  p(y | 6,p)  =  0.5  Zfn(l-20TC. .+0TQ. .0) 

ft  ~  ~1D  ~ 

-  (M2/2 )  £n2irp  -  ±  Z  ( (y  (i,  j ) -0Tz  (i ,  j  )  )  2  (2.15) 

ft  ~  ~ 

The  ML  estimates  0*  and  p*  are  obtained  by  maximizing  (2.15) 
with  respect  to  0  and  p  and  are  given  below: 


0*  -  Arg  min  {-0.5  Zfn(l-26TC. .  +  0TQ. .0) 

g  ft  ~~ij  ~~ij~ 

+  (N/2)fn  My(i,j)-0Tz(i,j)>2 


P*  =  -4  Z (y (i, j )-0*Tz  (i, j ) ) 2 

M  ft  ~  “ 


(2.16) 


(2.17) 


Since  the  log  likelihood  function  is  non-quadratic  in  0, 
the  estimation  involves  the  use  of  numerical  optimization 
methods  such  as  Newton- Rophson  [13]  which  are  computationally 
expensive.  We  given  an  iterative  method  which  yields  estimates 
close  to  the  ML  estimates  with  a  faster  convergence  rate.  The 
estimation  scheme  is  given  by 

Theorem  1:  Let  the  observations  y (i, j) , (i, j) 6Q  obey  the  R.F. 
model  in  (2.1)  characterized  by  6,p.  The  estimates  0", ~p  are 
given  by  the  following  iterative  scheme: 

0t+l  =  (Rt  -  J:f)“1(V  -  z^)  (2.18) 


and 


pt+l  “  "  ®t+l5(i,:’)) 

M  ' 


with 


-1 


®0  88  ?  Hi 


P0  ■  -^E(y(i,j)  -  QqZ  <± ,  3 ) )  2 
M  ft 


S  =  Ez(i,j)z  (i, j) 
ft 

U.=  Ey (i, j)z(i, j) 

~  ft 


V  =  0.5  (E( 


aij 


ij 


ij  i  j 

•ij  -  o.sui-2^.  ^9ijet] } 


and 


(2.19) 


(2.20) 

(2.21) 

(2.22) 

(2.23) 

(2.24) 


j  and  are  as  in  Lemma  1. 


The  proof  of  Theorem  1  is  given  in  the  appendix. 

Comments :  (1)  The  iterative  scheme  is  obtained  by  approximating 

the  determinant  term  up  to  quadratic  terms  using  a  logarithmic 
approximation . 

(2)  Without  losing  much  accuracy,  further  savings  in  compu¬ 
tation  can  be  achieved  by  letting  a^^  =  afc,  lsi,j&M 

where  a  =  {-4-  Ea. When  this  approximation  is  made,  the 
1  M  ft  1] 

vector  V  in  (2.22)  is  identically  equal  to  0  (due  to  a  property 

of  sums  of  exponentials)  and  the  following  computational  scheme 

results:  _  . 

S  1 

et+i  =  _(Rt“  ui>  (2.25) 

pt  pt  ~ 

pt+1  =  -^(y(i,j)  -  ft+15(i»j))2  (2.26) 

M 


where 


R  =  0 . 5  {  (~ 
at 


if  n  13 


EC. .cT. } 
at 


2 

2 


(2.27) 


and  S  is  as  in  (2.20). 

We  give  a  brief  discussion  regarding  the  decision  rule  for 
the  choice  of  appropriate  neighbor  sets.  Suppose  we  had  three 
sets  ,  n2,  and  of  neighbors  containing  m^,  m2,  and  n»3 


neighbors,  respectively.  Corresponding  to  each  N^,  we  write  the 
RF  model  as 


y(i,j)  =  E  0  k  0  y  [  (i,  j)+(k,£)  ]  +  /p“w(i,j)  (2.28) 

(Jc^)CNq  q'K'4  q 

0q,k,f*°  <k'*)tNq'  q  =  1,2,3 

Then  [8,10,11]  the  decision  rule  for  the  choice  of  appropriate 
neighbors  is:  choose  the  neighbor  set  if 

k*  =  arg  mintC.  >  where 
k  K 
T 

C.  =  {-E£n (1-20C.  •  +6TQ.  .0)  +  M2£np  +  M.£n  M2}  (2.29) 

K  q  - **  K 

Suppose  that  we  wish  to  also  include  unilateral  or  causal 
RF  models;  then  the  decision  statistic  reduces  to 

Ck  =  M2^np  +  M^n  M2  (2.30) 

This  follows  from  the  fact  that  the  Jacobian  of  the  transforma¬ 
tion  matrix  B(0)  from  noisy  variates  to  observations  is  unity 
[9,11].  Hence  the  model  selection  procedure  consists  of  comput¬ 
ing  Ck  for  different  models,  depending  on  whether  the  under¬ 
lying  models  are  causal  or  noncausal,  and  choosing  the  one  corre¬ 
sponding  to  the  lowest  G^. 


3.  Experimental  results 

We  describe  the  results  of  some  experiments  regarding  the 
generation  of  two-dimensional  patterns  and  estimation  schemes 
developed  in  the  previous  section. 

Experiment  1:  Synthetic  generation  of  two-dimensional  patterns. 
From  (2.4)  we  have,  for  the  observation  set  y  obeying  an 


RF  model. 


B(0)y  =  /p  w 


(3.1) 


The  synthetic  generation  is  then  done  by  assigning  some  arbi¬ 
trary  values  in  the  stationary  region  to  6  and  p  and  using  a 
pseudorandom  number  generator  to  form  the  vector  ui.  Since  the 
matrix  B(0)  has  a  block  circulant  structure,  Fourier  computa¬ 
tions  can  be  used  for  solving  y  in  (3.1)  [10].  Before  proceed¬ 
ing  further,  we  need  to  define  the  following  quantities:  denote 

2  .  . 
the  M  Fourier  vectors  fi;.,  by 

fi;j  =  column  (t  j  ,  Xitj, . . .  ,XiM-1t..) 

2  M™*  1 

tj  =  column (1,X ,... ,Xj  ±) ,  M-vector 
X±  =  exp[v^T  2ir(i-l)/M] 

The  synthetic  generation  scheme  then  is  as  follows: 


y  =  ./y.  .)  +  ai 


-13-13  13 


(3.2) 


where 


/p  r 

X.  .  =  — 5r  E  f.  .  U> 

1:  -1:  - 


u,.  «  (l-e'f,.),  i*i,j»M 


He  generate  the  vector  w  from  pseudorandom  numbers , 


generate  its  Fourier  sequence  {x^}  by  a  two-dimensional  FFT, 
and  finally  use  (3.2).  16  such  64x64  images  were  generated 

using  RF  models  with  different  neighbor  sets  and  parameters. 

The  gray  scale  values  of  the  images  were  corrected  to  lie  in 
the  range  0-63.  The  details  of  the  models  are  given  in  Table  1 
and  the  corresponding  images  are  shown  in  Fig.  1.  It  can  be 
seen  that  the  patterns  generated  are  quite  varied  and  some  of 
them  look  similar  to  real  textures.  Contrary  to  the  existing  [16] 
belief  that  autoregressive  RF  models  are  incapable  of  exhibit¬ 
ing  the  local  pattern  replication  attribute  considered  an  essen¬ 
tial  ingredient  of  texture,  some  of  the  windows  do  exhibit 
repetitive  patterns. 

We  use  matrix  notation  in  referring  to  windows  of  images. 

The  (1,1)  image  illustrates  the  idea  that  causal  models  are 
also  capable  of  accounting  for  some  periodic  patterns.  The 
parabolic  neighbor  set  seems  to  induce  vertically  oriented 
patterns  in  the  (1,2)  image.  The  (1,3)  and  (1,4)  windows  have 
identical  parameter  values  but  correspond  to  neighbor  sets  that 
are  related  (nearly  a  mirror  reflection)  and  result  in  similar 
diagonal  patterns  but  oriented  differently. 

Diagonal  neighbors  seem  to  induce  diagonal  patterns  as  in 
the  (2,2),  (3,4),  (4,1),  (4,3),  and  (4,4)  windows.  The  windows 
in  Fig.  1  are  only  typical  examples  and  more  interesting  patterns 
could  be  generated  by  varying  the  neighbor  sets  and  parameters. 


Experiment  2 :  The  role  of  parameter  values  in  the  structure 
of  patterns. 

To  illustrate  the  role  played  by  the  coefficients  in 
generating  the  two-dimensional  patterns,  we  consider  the  pattern 
corresponding  to  the  RF  model  N  =  { (-1,1) , (1,1) , (1,-1) ,  (-1,-1)  } , 
a  =  30.000,  p  =  1.1111  and  0_1  i  =  -i  =  and  i  =  _i 

a  .28.  The  values  of  the  parameters  tried  are  given  in  Table  2 
and  the  corresponding  pictures  in  Fig.  2.  Note  that  as  the 
parameters  are  varied  the  basic  pattern  is  still  retained  but 
the  "busyness"  of  the  pattern  is  varied.  Also  note  that  by 
changing  the  sign  of  the  diagonal  weights,  diagonal  patterns 
of  opposite  orientation  are  produced.  All  the  patterns  considered 
thus  far  were  generated  using  the  same  pseudorandom  number 
generator.  The  role  played  by  using  different  sets  of  pseudo¬ 
random  numbers  is  illustrated  later. 

Experiment  3:  To  test  the  usefulness  of  the  estimation  scheme 
and  the  choice  of  appropriate  neighbor  sets,  experiments  were 
done  with  two  synthetic  patterns .  The  true  model  corresponding 
to  the  first  pattern  is  defined  as  follows:  the  values  of  a  and 
p  are  30  and  1.1111  respectively,  the  neighbor  set  N  =  {(-1,0), 
(-1,1), (1,1), (1,0), (1,-1), (-1,-1)},  and  e_lfQ  =  eifQ  =  .12, 

0_2^2  *0^  *  .28,  0^  ^  =  0_i  =  -.14.  Using  the  model,  the 

synthetic  image  (1,1)  in  Fig.  3  was  generated.  But  for  correct 
inference  purposes  regarding  the  estimation  schemes,  the  original 


window  (values  not  scaled  for  display  purposes)  was  used. 

For  estimation  of  the  parameters,  the  sample  mean  of  the 
window  was  subtracted  and  the  iterative  scheme  developed  in 
(2.25-2.27)  was  used  for  different  RF  models.  The  test  statis¬ 
tic  defined  in  (2.29)  or  (2.30)  was  also  computed.  The 
actual  values  of  the  estimates  corresponding  to  different 
neighbor  sets  are  given  in  Table  3  and  the  corresponding  re¬ 
constructed  images  are  given  in  Fig.  3. 

Table  3  shows  that  the  estimated  values  corresponding  to 
the  true  neighbor  set  are  close  to  the  true  values.  The  so- 
called  least  square  estimates  corresponding  to  the  true  neigh¬ 
bor  set  are  inefficient  compared  to  the  estimates  obtained  by 
the  scheme  developed  here.  Note  that  when  extra  neighbors 
are  added,  the  corresponding  parameter  values  are  very  small. 

The  decision  statistics  corresponding  to  the  models  considered 
are  tabulated  in  Table  4  and  the  decision  rule  picks  up  the 
true  model. 

Fig.  3  shows  the  pictures  corresponding  to  the  different 
neighbor  sets  using  an  identical  array  of  noise  variables  in 
Table  3.  The  (1,2)  window  corresponding  to  the  least  square 
estimate  of  the  true  model  is  a  poor  reproduction  of  the  original 
in  (1,1).  On  the  other  hand,  the  (1,3)  window  corresponding 
to  the  approximate  ML  estimates  developed  here  is  very  close 


to  the  original.  The  reproductions  corresponding  to  neighbor 
sets  {(-1,0),  (1,0),  (0,-1)}  and  { (-1,0) ,  (0 ,1) ,  (1,0)  ,  (0,-1)}  are 
very  poor.  The  windows  (2,4),  (3,1),  and  (3,2)  look  close  to 
the  originals  as  the  corresponding  models  include  the  original 
neighbor  set.  But  the  decision  rule  suggested  here  correctly 
eliminates  these  models. 

Since  in  real  world  applications  the  procedure  of  using 
identical  noise  variables  for  the  reconstructed  pictures  as  for 
the  originals  is  unrealistic,  the  natural  question  is  how  sensi¬ 
tive  are  these  patterns  to  different  sequences  of  noise  varia¬ 
bles.  To  answer  this  question  synthetic  generation  of  patterns 
were  done  using  different  random  number  sequences  and  the  true 
model  as  in  Expt.  3,  but  with  estimated  parameters.  The  result¬ 
ing  patterns  are  shown  in  Fig.  4.  Note  that  the  variations  in 
the  patterns  are  noticeable  if  one  taks  a  close  look,  but  the 
basic  patterns  are  retained  in  all  the  windows.  (Window  (1,3) 
corresponds  to  the  random  numbers  used  in  Fig.  3  and  consequently 
is  identical  to  window  (1,3)  of  Fig.  3.) 

Table  5  and  Fig.  5  correspond  to  experiments  with  the  syn¬ 
thetic  pattern  generated  by  the  causal  model  {  (-1,0) ,  (0,-1) , 
(-1,-1)}  with  a  *  30.8870,  p  =  0.1087,  0^  Q  =  .9704,  eQ  _1  = 
.9735  and  8_^  =  -.9686.  The  estimates  corresponding  to 

different  RF  models  are  given  in  Table  5.  While  estimating  the 
parameters  for  causal  models  the  least  square  estimates  them¬ 
selves  were  treated  as  the  approximate  ML  estimates,  though  the 


underlying  patterns  were  produced  using  a  torus  structure.  (The 
error  in  approximation  due  to  the  torus  representation  used 
for  the  causal  models  is  negligible.)  Note  that  the  decision 
rule  (see  Table  6  for  the  test  statistics)  picks  up  correctly 
the  causal  model  compared  to  other  noncausal  models.  This  shows 
that  it  is  not  true  that  noncausal  models  are  always  superior  to 
causal  models. 

Fig.  5  shows  the  windows  constructed  using  the  models  in 
Table  5.  The  quality  of  reproduction  using  the  causal  model 
is  markedly  superior  to  that  using  the  semicausal  model 
{  (1,4) , (2 ,1) }  or  noncausal  model  (2,2).  (In  fact,  the  noncausal 
model  t  (-1,0) , (0,1) , (1,0) , (0,-1) }  yielded  estimates  in  the 
nonstationary  range.)  Since  the  extra  neighbors  used  in  the 
models  corresponding  to  (1,4)  and  (2,1)  have  very  small  values 
these  patterns  look  very  similar.  As  more  and  more  noncausal 
members  are  added  (windows  (3,1)  and  (3,2))  the  reproduction  is 
somewhat  better.  Windows  (2,3)  and  (2,4)  correspond  to  causal 
neighbor  sets  that  include  the  original  model  and  hence  the 
quality  of  reconstruction  is  good. 

The  variations  produced  in  the  pattern  corresponding  to 
different  random  number  sequences  (considered  in  Fig.  4)  are 
shown  in  Fig.  6.  Note  that  the  basic  patterns  are  retained. 


i 


4.  Discussion 

We  have  considered  the  problem  of  estimating  the  unknown 
parameters  of  an  RF  model  and  the  choice  of  appropriate  neigh¬ 
bors.  The  iterative  estimation  scheme  yields  estimates  that 
are  close  to  the  ML  estimates.  The  problem  of  statistical 
inference  of  stationary  RF  models  has  been  previously  considered 
by  Whittle  [12]  and  Larimore  [9] .  Since  in  [12]  representation 
on  a  square  lattice  was  used,  the  evaluation  of  the  determinant 
of  the  transformation  matrix  is  difficult,  and  approximate 
methods  using  power  series  expansion  of  the  spectral  density 
were  used.  The  computation  of  the  determinant  can  be  conveni¬ 
ently  done  in  the  transform  domain  as  in  [9],  but  the  method 
involves  the  assumption  that  the  Fourier  transforms  of  y ( • ) 
are  uncorrelated,  which  is  true  only  for  an  infinite  image. 

By  using  the  torus  representation  for  finite  images,  as  done  in 
this  paper,  the  determinant  can  be  explicitly  evaluated  and  the 
exact  likelihood  function  can  be  written  down.  This  leads  to 
computations  in  the  spatial  domain  as  against  the  transform 
domain. 

The  iterative  scheme  suggested  here  yields  estimates  that 
are  close  to  the  true  parameters,  with  less  computational  effort 
compared  to  numerical  optimization  methods  such  as  Newton-Raphson, 


etc. 


We  have  also  implemented  a  decision  rule  for  the  choice  of 
neighbors  which  correctly  picks  up  the  true  model  against  many 
competing  models. 

To  illustrate  that  RF  models  are  indeed  capable  of  generat¬ 
ing  a  wide  variety  of  patterns ,  examples  of  synthetic  genera¬ 
tion  results  have  been  presented.  Most  of  the  patterns  possess 
the  local  pattern  replication  property  which  is  considered  to 
be  an  essential  ingredient  of  textures. 

We  have  illustrated  the  theory  using  synthetic  patterns. 
Currently,  work  is  under  progress  in  testing  the  estimation 
scheme  with  real  textures. 


References 


1.  A.  Habibi,  "Two-dimensional  Bayesian  estimate  of  images", 
Proc.  IEEE,  Vol.  60,  pp.  878-883,  July  1972. 

2.  D.  P.  Panda  and  A.  C.  Kak,  "Recursive  least  squares  smooth¬ 
ing  of  noise  in  images",  IEEE  Trans,  on  Acoustics,  Speech 
and  Signal  Processing,  Vol.  ASSP-25,  pp.  520-524,  Dec.  1977. 

3.  A.  K.  Jain  and  J.  R.  Jain,  "Partial  difference  equations 
and  finite  difference  methods  in  image  processing — part  2 : 
image  restoration",  IEEE  Trans,  on  Automatic  Control, 

Vol.  AC-23,  pp.  817-833,  Oct.  1978. 

4.  E.  J.  Delp,  R.  L.  Kashyap  and  0.  R.  Mitchell , "Image  data 
compression  using  autoregressive  time  series  models". 

Pattern  Recognition,  Vol.  11,  pp.  313-323,  Dec.  1979. 

5.  A.  K.  Jain  and  E.  Angel,  "Image  restoration,  modeling  and 
reduction  of  dimensionality",  IEEE  Trans,  on  Computers, 

Vol.  C-23 ,  pp.  470-476,  May  1976. 

6.  K.  Degucni  and  I.  Morishita,  "Image  coding  and  reconstruc¬ 
tion  by  two-dimensional  optimal  linear  estimation",  Proc. 

4th  IJCPR,  pp.  530-532,  Nov.  1978. 

7.  J.  T.  Tou,  Pictorial  feature  extraction  and  recognition 
via  image  modeling,  Computer  Graphics  Image  Processing, 

Vol.  12,  pp.  376-406,  April  1980. 

8.  R.  L.  Kashyap,  "Univariate  and  multivariate  random  field 
models  for  images”.  Computer  Graphics  Image  Processing, 

Vol.  12,  pp.  257-370,  March  1980. 

9.  W.  E.  Larimore,  "Statistical  inference  on  stationary  random 
fields",  Proc.  IEEE,  Vol.  65,  pp.  961-970,  June  1977. 

10.  R.  L.  Kashyap,  "Random  field  models  on  torus  lattices  for 
finite  images"  (submitted) . 

11.  R.  L.  Kashyap,  R.  Chellappa,  and  N.  Ahuja,  "Decision  rules 
for  choice  of  neighbors  in  random  field  models  of  images" 

(to  appear  in  Computer  Graphics  Image  Processing) . 

12.  P.  Whittle,  "On  stationary  processes  in  the  plane", 
Biometrika,  Vol.  41,  pp.  434-449,  1954. 


r 


13.  N.  K.  Gupta  and  R.  K.  Mehra,  "Computational  aspects  of 
maximum  likelihood  estimation  and  reduction  in  sensitivity 
function  calculations",  IEEE  Trans,  on  Automatic  Control, 
VQl.  AC- 19,  pp.  774-783,  Dec.  1974. 

14.  R.  L.  Kashyap,  "A  Bayesian  comparison  of  different  classes 
of  dynamic  models  using  empirical  data",  IEEE  Trans .  on 
Automatic  Control,  Vol.  AC-22,  pp.  715-727,  Oct.  1977. 

15.  R.  L.  Kashyap  and  R.  Chellappa,  "Estimation  methods  for 
random  field  models  of  images"  (in  preparation) . 

16.  J.  W.  Modestino,  "Texture  discrimination  based  upon  an 
assumed  stochastic  texture  model",  TR-79-3,  Electrical  and 
Systems  Engineering  Dept.,  Rensselaer  Polytechnic  Institute, 
Troy,  New  York,  July  1979. 


Proof  of  Lenuna  1: 


Appendix 


From  (2.11), 

In  det  B  (0)  =  Zln  (l-0T'i' .  . )  (1) 

Using  the  definition  of  4* .  .  . 

(l-eV.)  =  (i-0Tci j  -  (2) 

where 

cjj  =  {cos  A0l(i-l)k+(j-l)flf  (k,*)fcN} 

and 

=  (sin  XQ  [  (i-l)k+(j-l)£] ,  (k,£)  fcN} 

From  (2) , 

_*TS 

{n(l-eV)  =  J  £n[  (1-0TC.  . ) 2  +  (0TS.  .)2]  +  tan~1(---;^) 

J  ~  ■vl-)  ~  ~13  i-e1cij 

Since  Zln  (1-0T4I .  . )  is  real 
SI  ~ 

In  det  B (0)  =  E-£n  (l-0T'i’ .  .) 

si  ~ 

=  O.5[E£n(l+0T(C.  .Crf  .  +  S .  .  S*  ) -26TC .  . )  (4) 

~  ~i3~ij  -13-13  -  -13 


Defining 


Q. .  =  C. .C.?  +  S. .S.T 
~13  -1]-13 


(3) 


proves  Lemma  1. 


Proof  of  Theorem  1:  Prior  to  proving  Theorem  1  we  need  the 
following  lemma  which  is  proved  subsequently. 

Lemma  2 : 


Lena-e^j)  =  K(aij)-yTe  +  0tr0  (5) 

where  the  vector  V  and  the  matrix  R  are  as  in  (2.22)  and  (2.23) 
and 

K (a± j )  =  0.5{Un  a± (1^ . ) /a±  .  -  (l-ai  ) /2a* ^  } 

From  (2.15)  and  Lemma  2, 

£np(y|0,p)  =  K (a± j )  -  VT0  +  0TR0  - (M2/2) £n2np 

-  A.  E (y (i, j )  -  0Tz(i,j))2  (6) 

zp 

Differentiating  w.r.t.  0  and  p  and  equating  to  zero. 


and 


-V  +  R0  +  \  Ez (i, j) (y (i, j)-zT(i, j) 0) 

p  £2~ 


=  0 


(7) 


P  =  -K  £(y(i,j)  -  0z  (i,  j) ) 2 

m  n 


(8) 


Solving  (7),  and  defining 

S  =  Ez (i, j)zT(i, j)  ,  (9) 

(T 

U  =  Ez(i, j)y (i, j)  (10) 

~ x 

we  have 

0  =  (H  -  z  S)_1(V  -  i  Ux)  (11) 

P  -  P 

From  (11)  and  (13)  we  obtain  the  following  iterative  scheme: 


i  , 

'i  * 


Defining 


V 

and 


0.51 - 
ft  aij 


R  =  0.5S[  -  4 

ft  aij  a‘ 

we  arrive  at  Lemma  2. 


NUMBER  NEIGHBOR  SET  d  p  COEFFICIENTS 


Table  1  (cont'd 


NUMBER 


(1,1) 


(1,2) 


(1,3) 


(1,4) 


(2,1) 


(2,2) 


(2,3) 


(2,4) 


(3,1) 


(3,2) 


(3,3) 


(3,4) 


(4,1) 


(4,2) 


(4,3) 


(4,4) 


el,l  =  °-l,-l 

ei.-i  =  e-i,i 

.28 

-.14 

.28 

-.10 

.28 

-.06 

.30 

-.10 

.32 

-.10 

.34 

-.10 

.36 

-.10 

.38 

-.10 

-.14 

.28 

-.10 

.28 

-.06 

.28 

-.10 

.30 

o 

H 

• 

1 

.32 

O 

H 

• 

1 

.34 

O 

H 

• 

1 

.36 

O 

H 

• 

1 

.38 

\i 


i 

>0 

rH 

p 

• 

c 

4) 

N 

tr  4i 

0 

'O 

Q 

n  p 

Oi 

o 

o 

2 

n  o 

e 

as 

P  6 

4)  P 

«  -H 

p 

4) 

w 

<0  P 

P  P  3 

P 

<u  u> 

o  c 

P 

« 

£h 

P  0) 

O  -H 

P 

(0 

VO 

H 

CM 

z 

** 

o 

in 

O 

w 

rH 

00 

i" 

rH 

CO 

M 

CN 

• 

CN 

•<r 

• 

VO 

o 

rl 

| 

• 

rH 

» 

cn 

H 

• 

II 

II 

• 

II 

• 

h 

II 

rH 

rH 

II 

rH 

II 

IP 

O 

1 

•k 

o 

rH 

W 

* 

rH 

* 

| 

o 

H 

H 

i 

rH 

rH 

% 

o 

CD 

1 

<D 

CD 

1 

o 

II 

CD 

II 

ii 

CD 

CD 

O 

II 

rH 

o 

II 

II 

* 

rH 

1 

w 

* 

rH 

rH 

rH 

«k 

rH 

1 

H 

rH 

1 

rH 

O 

CD 

<D 

CD 

CD 

CD 

CD 

H  O  «H 
I  *  I 


■  v  -/  »  ^  t  ■  >«>pT5! 


NUMBER  NEIGHBOR  SET  a  p  COEFFICIENTS 


NUMBER  I  NEIGHBOR  SET  [  a  I  p  I  COEFFICIENTS 


NUMBER 


TEST  STATISTIC 


(1.2) 

-8900.0 

(1.3) 

3883.16 

(1,4) 

-8663.66 

(2,1) 

-8657.67 

(2,2) 

(non- stationary ) 

(2,3) 

-8892.0 

(2,4) 

| 

-8884.6 

(3,1) 

-6102.58 

(3,2) 

-6866.006 

Figure  1.  Examples  of  synthetic  generation  of  images 
using  RF  models  (see  Table  1) . 


Figure  2 


Patterns  produced  by  models  with  the  same  neighbor 
set  N={  (-1,1)  ,  (1,1)  ,  (1,-1)  (-1,-1)  },  a=30 . 034  ,  p=1.1240 
but  with  different  sets  of  values  for  the  coefficients 
(see  Table  2) . 


Figure  3.  Reconstruction  of  images  corresponding  to  the  true 
model  N={  (-1,0) ,  (-1,1) , (1,1) , (1,0) , (1,-1) ,  (-1,-1) }. 
See  Table  3  for  details. 


Figure  5.  Reconstruction  of  images  corresponding  to  the  true 
model  N={ (-1,0) , (0,-1) , (-1,-1) } .  See  Table  5  for 
details. 


Figure  6. 


Patterns  produced  using  different  sequences  of  random 
numbers  and  the  model  N={ (-1,0) , (0,-1)  (-1,-1) } , 


a=31. 096 ,  p=0 . 1136 ,  8  ,  Q=.9772, 
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